#include "HadronicWReconstructor.h"
#include <gsl/gsl_combination.h>
HadronicWReconstructor::HadronicWReconstructor()
{

}

HadronicWReconstructor::HadronicWReconstructor ( EvtData* e,int s)
{
scheme=s;
evtdata=e;
WContainer = new std::vector<WCandidate>;

}

HadronicWReconstructor::~HadronicWReconstructor()
{
WContainer->clear();
delete WContainer;
}

WCandidate HadronicWReconstructor::getW()
{
 W=WContainer->at(0);
return W;
}

bool HadronicWReconstructor::MinimizebyDeltam ( WCandidate w1, WCandidate w2)
{
return ( fabs(w1.theW.P4().M()-MW )<fabs(w2.theW.P4().M()-MW ));
}
bool HadronicWReconstructor::MaximizebyPt ( WCandidate w1, WCandidate w2)
{
return ( (w1.j1.P4()+w1.j2.P4()).Pt() > (w2.j1.P4()+w2.j2.P4()).Pt() );
}


unsigned int HadronicWReconstructor::getNWcand()
{
return NWcand;
}

void HadronicWReconstructor::Run()
{
  //std::cout<<"HadronicWReconstructor::jets_size="<<evtdata->getJets()->size()<<std::endl;

  WContainer->clear();
  if(scheme==0 || scheme==1){
  WCandidate Wboson;
   //Select any 2 jets from the remaining and loop over different 2 jets combinations 
  gsl_combination * comb;
  comb = gsl_combination_alloc(evtdata->getJets()->size(),2);
  gsl_combination_init_first(comb);
  //recoEvt DquarkPair;
 
 while (true){	//Jet combination loop
	 //get indices of W constituent jets
	 unsigned int i1=gsl_combination_get(comb,0);
	 unsigned int i2=gsl_combination_get(comb,1);
	 particle Wconstituent1(evtdata->getJets()->at(i1));
	 particle Wconstituent2(evtdata->getJets()->at(i2));
	 Wboson.j1=Wconstituent1;
	 Wboson.j2=Wconstituent2;
         Wboson.i1=i1;
         Wboson.i2=i2;
         particle Wparticle;
         Wparticle.setP4((Wconstituent1.P4()+Wconstituent2.P4()));
         Wboson.theW=Wparticle;
	 double delta_m=fabs((Wboson.j1.P4()+Wboson.j2.P4()).M()-MW);
	 Wboson.deltaM=delta_m;
         if( Wparticle.P4().M()>60 && Wparticle.P4().M()<100 && Wconstituent1.P4().DeltaR(Wconstituent2.P4())<1.75){
        WContainer->push_back(Wboson);
         }
     if(gsl_combination_next(comb)!=GSL_SUCCESS) break;
 }//endof combination loop
gsl_combination_free(comb);
NWcand=WContainer->size();
//never forget to remove these jets !!!
}
if(scheme==0){
std::sort(WContainer->begin(),WContainer->end(),MinimizebyDeltam);
}
else if(scheme==1){
std::sort(WContainer->begin(),WContainer->end(),MaximizebyPt);  
}
//single W jet case
/*
else if(scheme==2){
 WCandidate Wboson;
 particle Wparticle;
 for(unsigned int i=0;i<evtdata->getJets()->size();i++){
   double mjet=evtdata->getJets()->at(i).P4().M();
   if(mjet<60 || mjet>100)continue;
     Wparticle.setP4(evtdata->getJets()->at(i).P4());
     Wboson.theW=Wparticle;
     Wboson.i1=i;
     Wboson.i2=-1;
     Wboson.j1=Wparticle;
     Wboson.j2=Wparticle;
     //std::cout<<Wboson.theW.P4().Pt()<<std::endl;
     WContainer->push_back(Wboson);
 }//endof jet loop 
NWcand=WContainer->size();
std::sort(WContainer->begin(),WContainer->end(),MinimizebyDeltam);
}//endof scheme
*/
else{
  
  std::cout<<"ERROR***No such scheme "<<std::endl;
  exit(-8);
}
return;

}
